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A new computational method is presented for study suspensions of charged soft particles under¬ 
going fluctuating hydrodynamic and electrostatic interactions. The proposed model is appropriate 
for polymers, proteins and porous particles embedded in a continuum electrolyte. A self-consistent 
Langevin description of the particles is adopted in which hydrodynamic and electrostatic interac¬ 
tions are included through a Green’s function formalism. An Ewald-like split is adopted in order 
to satisfy arbitrary boundary conditions for the Stokeslet and Poisson Green functions, thereby 
providing a formalism that is applicable to any geometry and that can be extended to deformable 
objects. The convection-diffusion equation for the continuum ions is solved simultaneously consid¬ 
ering Nernst-Planck diffusion. The method can be applied to systems at equilibrium and far from 
equilibrium. Its applicability is demonstrated in the context of electrokinetic motion, where it is 
shown that the ionic clouds associated with individual particles can be severely altered by the flow 
and concentration, leading to intriguing cooperative effects. 


There is considerable interest in understanding the 
structure and dynamics of suspensions of charged parti¬ 
cles over multiple length scales, both at equilibrium and 
far from equilibrium. Examples include DNA or protein 
flow in mocrofluidic devices, in cellular environments, or 
colloidal self-assembly in external fields. Beyond any di¬ 
rect interaction (van der Waals or electrostatic) between 
particles, the motion of a particle in solution induces 
important hydrodynamic and electrostatic interactions 
that some times compete against each other, leading to 
electro-osmotic or electro-kinetic phenomena that remain 
poorly understood. 

Resolving the dynamics of solvents or charged species 
over short and long length scales remains a significant 
challenge. The central question is how to evolve these 
systems, while preserving molecular resolution of discrete 
macromolecules or colloids, adopting continuum descrip¬ 
tions for solvent and ions (see Figure [T]). Past attempts 
have primarily relied on Lattice Boltzmann (LB) [l|. 
Stochastic Rotational Dynamics (^RD) [I, Q, Stokesian 
Dynamcis (SD) Q, Ewald sums 0 and the General ge¬ 
ometry Ewald-like method (GgEm) Q for the hydro- 
dynamic evolution, i.e. the momentum equations. On 
the other hand, solution of charges have been treated 
by Ewald-based methods 0 or Poisson-Boltzmann solu¬ 
tions [1, 0. A notable exception, where both interac¬ 
tions are solved simultaneously, is the Smoothed Profile 
Method (SPM) Eliii. In SPM, discrete particles are 
included into the continuum formalism through smooth¬ 
ing functions (see the review of particle-based methods 
by T. Yamaguchi et. al [ill). LB and SRD methods 
rely on collision operators to evolve fluid dynamics, thus 
imposing theoretical limits at zero inertia or strict incom¬ 
pressibility (Ma = Re = 0). In addition, they present 
limitations regarding the size of the systems that can 
be studied. A salient limitation of SPM and LB is the 
mesh dependency that naturally develops at finite con¬ 
centrations of the discrete entities. The method that we 
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FIG. 1. Discrete charged soft-particles embedded in an elec¬ 
trolyte solvent. (A) A level of description that requires to 
treat soft-particles, ions and fluid as discrete entities. (B) 
The proposed method provides resolution for the interesting 
discrete entities whereas a continuum description for the elec¬ 
trolyte solvent. 


propose in this work, provides a way of resolving these 
problems. 

We are presenting a novel theoretical method that re¬ 
solves the dynamical coupling between discrete charged 
soft-particles and continuum electrolyte solvents. It is 
presented as a generalization of the GgEm approach, 
which can be used for bulk and confined systems, both 
at equilibrium and far from equilibrium. The dynam¬ 
ics of the soft-particles follows a Fokker-Planck equation 
for the probability density, resulting in a Stochastic evo¬ 
lution equation. The solvent and the continuum ions, 
on the other hand, are evolved according to momentum 
and mass and balances, including Nernst-Planck diffusion 
mechanisms (see Fig. [T]). 

We consider a collection of Np soft-particles, with 
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charge eZy and hydrodynamic radius a, suspended in 
a solvent that includes Nj charged species (continuum 
ions). The particles are “soft” in the sense that they are 
permeable to the continuum ions and the fluid flow. This 
model is suitable for coarse-grained models of polymers, 
proteins and other biological systems. Soft-particles and 
continuum ions contribute to a charge density, defined as 

Ni Np 

p(x) =F'^ZjCj{x.) + ^ Z,,e6{x - x^), (1) 

j = l U=1 

where F = eN^ is Faraday’s constant {Na is Avo- 
gadro’s number), zj is the valence of the continuum ions 
{j = 1,..., A/"/), Cj is the concentration of the continuum 
species, Zp is the soft-particle’s valence {u = l,...,A/'p) 
and e is the elementary charge. Where the soft-particles, 
at this point, are considered point-charges. If electroneu¬ 
trality is not satisfied at a local level, the charge density 
will drive an electrostatic potential given by the solution 
of Poisson’s equation: 

V^0(x) = -p(x)/eoe, (2) 

where cq is the vacuum permittivity and e is the solvent 
relative permittivity. The electric field (E(x) = — V0(x)) 
drives an electric force on the ions. Electric forces on the 
continuum ions and the total forces on the soft-particles 
define a force density given by 


Ni Np 

p(x) = F^2 ;jP(x)E(x) + ^f^(5(x-x^), (3) 

j = l U=1 

where represents the total non-Brownian and non¬ 
hydrodynamic force acting on particle z/. Neglecting in¬ 
ertia {Re = 0), the solvent velocity can be written as 
v(x) = vo(x) -h u(x), (where vo(x) is the unperturbed 
velocity and u(x) is the velocity perturbation), and it is 
given by the solution of a Stokes system of equations. 
The velocity perturbation is driven by the force density, 
i.e. 


- Vp(x) + 7^V^u(x) = -p(x), V • u(x) = 0, (4) 

where r] is the solvent viscosity. 

The evolution of the ions within the solvent follows a 
species balance with the total flux defined as a sum of 
convection and diffusion fluxes. The latter are given by 
the Nernst-Planck diffusion, resulting in: 

DjZj{e/kBT) ■ V<i>] , 

where Dj is the diffusion coefficient of ion j, ks is Boltz¬ 
mann’s constant, and T is the absolute temperature. 

Each of the Np discrete ions behaves, for the moment, 
as a point-force and a point-charge. The discrete ions 


have a hydrodynamic radius a, and interact via a re¬ 
pulsive excluded volume potential. Neglecting inertia 
(Reynolds, Re = 0), for each discrete ion, = 1, ...,Np, 
the force balance requires -T H- H- = 0, 

where, for bead z^, is the hydrodynamic force, is the 
particle-particle excluded volume force, is the electro¬ 
static force, is the Brownian force and are any other 
potential forces that may arise in the system IJ, Q. 
The dynamics of the discrete ions in the solvent is de¬ 
scribed by the probability density in configuration space. 
The diffusion equation for the configurational distribu¬ 
tion function has the form of a Eokker-Planck equation, 
which corresponds to the following system of stochastic 
differential equations of motion for the discrete ion posi¬ 
tions: 


dR = 


Vo + 


kpT 


D F 


A 


D 


dt + V2B-dW. (6) 


Here, R is a vector containing the 3Np coordinates 
of the soft-particles, and where R^y = x^y denotes the 
Cartesian coordinates of particle z/. The vector Vq, of 
length 3A/’p, represents the unperturbed velocity field, 
with Vo,iy = vo(xjy). The vector F has length 3A/’p, 
with Fjy = fjy being the total non-Brownian, non¬ 
hydrodynamic force acting on bead u. Finally, the in¬ 
dependent components of dW are obtained from a real¬ 
valued Gaussian distribution function with zero mean 
and variance dt. The diffusion tensor D (or mobility 
tensor) is a 3Np x 3Np tensor. It may be separated 
into the Stokes drag and the hydrodynamic interaction 
tensor, = kpT [SSj,^ + (1 - where 

^ is a 3 X 3 identity matrix and is the Kronecker 
delta. The Brownian perturbation, dW, is coupled to 
the hydrodynamic interactions through the fluctuation- 
dissipation theorem: D = B • B^. 

The characteristic variables for the system are set by 
the soft-particle: hydrodynamic radius, a, for length, 
particle diffusion time, C,o? jkpT^ for time (where C = 
GTTTja is the drag coefficient), e/Aneoea for the electro¬ 
static potential, and the elementary charge e for the 
charge. There are two scales for velocity: one for the un¬ 
perturbed velocity field vq and one for the velocity fluctu¬ 
ations Uc = kpTjC^a. The uniform concentration for one 
of the species, Co, is used as the characteristic concentra¬ 
tion of ions within the solvent. Therefore, Pe^- = voa/Dj 
defines a Peclet number for species j based on the im¬ 
posed flow field Vo and = (^Dj/kpT is the ratio be¬ 
tween particle and continuum-ion diffusion coefficients. 
The ratio between electrostatic forces and thermal forces 
defines the so-called Bjerrum length, Xp = /47reoe/cpT, 
and the ionic strength, I = ^ ’ defines the so- 

called Debye length, = 2NAe^I/eoekpT. 

The key feature of GgEm-like methods is to decompose 
the charge-density and force-density expressions into a lo¬ 
cal (free-space) contribution and a global (bounded) con¬ 
tribution, in analogy to Ewald-like methods. The “local” 
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densities are defined by 

Np 

^ 2 ;^. [6{x - Xj,) - gsix - x^)], 

( 7 ) 

y] e [<^(x - Xi.) - gH{x - Xi.)], 

U=1 

which produces a local contribution to the electrostatic 
potential 0/(x) and the velocity perturbation u/(x). The 
“global” densities are given by 


Pi(x) = 
Pj(x) = 


Ni Np 

j = l V=1 

Ni Np 

Pg{x) = Aft y] ZjpjCj'E + y] (x - x^)], 

j=l V=1 

and are responsible for the global contribution of the 
potential 0^(x) and velocity perturbation u^(x). The 
linearity of the Poisson and Stokes equations implies 
that 0(x) = 0z(x) + and u(x) = u^(x) + 

u^(x). The screening functions, gE,H{'^)^ satisfy 
/all space = 1. The local contributions, 

0/(x) and u/(x), are calculated analytically assuming 
an unbounded domain, only considering neighbor par¬ 
ticles. For Poisson’s equation, the screening func¬ 
tion is a Gaussian: ^£;(r) = {a^^ \ For 
the Stokes equations it is a modified Gaussian 

[5/2 — a‘^r‘^] (The local cal¬ 
culation scheme is described in the SI). The global con¬ 
tributions are found numerically, requiring that (j)i + (^g 
and Ui^Ug satisfy appropriate boundary conditions. For 
example, homogeneous Dirichlet boundary conditions re¬ 
quire (pgi'x.) = — 0/(x) and u^(x) = — u^(x). For prob¬ 
lems with periodic boundary conditions, Fourier tech¬ 
niques are used to guarantee the periodicity of the global 
contributions. The periodicity for local contributions is 
obtained through a minimum image convention. A nu¬ 
merical implementation for periodic domains is described 
in the SI. 

The point-source regularization (i.e. soft-particle) is 
implemented with the same screening functions from 
GgEm introducing two additional length scales and 
. These length scales are related to the hydrodynamic 
radius a. For the point-force (hydrodynamics), this is 
achieved by limiting the maximum velocity on the fluid 
driven by the regularized-force, which at i?e = 0 is given 
by Stokes’ law. The regularization for the point-charge 
is achieved by distributing the total charge throughout 
the particle, thereby “confining” the charge to the parti¬ 
cle volume. The regularization parameters are = 3/a 
and = y^/3a. 

The method is validated by relying on several approxi¬ 
mate analytical solutions for soft-particles. Most of these 


FIG. 2. Ionic cloud dimensions, in the neutral (Ln) and 
applied field direction (Le), as a function of the applied 
field for different Debye lengths (salt concentration). The 
soft-particles are embedded in a charged solvent and the 
applied field deforms the continuum ionic cloud. Blue iso¬ 
concentration surfaces depict the counter-ions surrounding 
the negatively charged soft-particles. 


solutions correspond to equilibrium conditions, and re¬ 
sult from the linearized Poisson-Boltzmann approxima¬ 
tion (i.e. low potentials). Solutions for the electrostatic 
potential and the electrostatic energy between particles 
are included in the SI. For instance. Fig. [2] illustrates how 
ionic clouds that surround the soft-particles are deformed 
when an external electric field is applied, as a function of 
the salt concentration. 

Non-equilibrium, finite concentrations, interacting 
ionic clouds and fluctuating hydrodynamic interactions 
are built in the NP-GgEm. To illustrate the behavior of 
these types of systems. Brownian soft-particles are simu¬ 
lated in a solvent at different (salt concentrations). 
Diffusing soft-particles of charge Zj^ = —1 and Zj^ = —2 
are simulated in a periodic box of size L = 20a. Note that 
global electroneutrality is enforced through the contin¬ 
uum ions. Soft-particle concentration is defined in terms 
of an effective volume fraction: (j)E = 4:Np7ra^/3L^. Fig¬ 
ure [3] provides representative results for NP-GgEm. In 
the figure, average counter- and co-ion concentration and 
RMS velocity profiles are shown as a function of distance 
from the soft-particle, for Xd = 2a. At equilibrium, the 
2D concentration profiles suggest that they follow PB- 
like behavior. However, instantaneous concentration iso¬ 
surfaces (ion clouds) demonstrate how, even at low con¬ 
centrations, the ion clouds merge and interact with each 
other, being deformed by the Brownian motion of the 
particles. The figure also includes the RMS velocity per¬ 
turbation (normalized by the number of particles) as a 
function of the distance from the soft-particle at to dif- 
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FIG. 3. (A) Average concentration and velocity fluctuations 
as a function of the distance from the soft-particle, r/a. Red 
and blue lines are for co- and counter-ion concentration pro¬ 
files, respectively; black lines depict RMS velocity profiles. 
(B) Instantaneous counter-ion clouds (blue) and (C) velocity 
streamlines (green). The soft-particles have negative changed 
and they are embedded in a solvent with Az:) = 2a at an ef¬ 
fective volume fraction 0e = 5.24 x 10“^ {N = 10). 


FIG. 4. (A) Diffusion coefficient as a function of the salt 

concentration (Debye length, Xd), soft-particles charge, z/, 
and effective volume fraction, 0^;. The diffusion coefficient 
is normalized with the particle diffusion coefficient Da = 
ksTlGiTija. (B) Electrophoretic mobility, /i, as a function 
of the effective volume fraction and salt concentration for a 
weak (blue lines) and a strong (black lines) imposed fields, 
Eq. (G) Electrophoretic mobility at zeroth concentration, /xo, 
as a function of the Debye length for the two applied fields in 
(B). The blue line represent Ohshima’s approximate solution 
for non-fluid-penetrated soft-particles [Iq . 


ferent concentrations. The flow streamlines (green) are 
also plotted. The RMS velocity exhibits a slow decaying 
tail, accounting for the long-range behavior of HI, and it 
increases non-linearly with concentration. 

One key property of interest in charged colloidal sys¬ 
tems is the diffusion coefficient, and its dependence on 
the concentration, charge and salt. To illustrate the 
behavior of this property, and the usefulness of the 
approach proposed here, in Fig. |4] we show the soft- 
particle diffusion coefficient as a function of the Debye 
length for two different concentrations and charges. The 
simulations were evolved from 200 to 1000 characteris¬ 
tic particle-diffusion times, thereby providing sufficient 
statistics. In general, the diffusion coefficient decreases 
by increasing soft-particle concentration and charge. For 
high salt (low Xd)^ the diffusion coefficient approaches 
the infinite dilute and non-charge system. As the salt 
concentration increases, electrolyte drag and collective 
motion drive a non-monotonic behavior, decreasing mo¬ 
bility at first and then increasing it. This behavior of 
diffusing charge particies in a charge soivent was aiso 
been observed in other systems |15l4l7j. In the figure, 
the eiectrophoretic mobility, jj. = u/Eq^ for soft-particles 
of charge = — 1 is also shown. The figure inciudes 
two appiied eiectric fieids: a weak Eq = 0.1 Ec and 
a strong Eq = 100Ec (where the characteristic fieid 
is Ec = e/47reoea^). We recaii that fluid can pene¬ 
trate the soft-particies, causing the eiectrophoretic mo¬ 
bility to change with the applied electric field. The 
electrophoretic mobility at zero concentration, /io, and 
an approximate solution for non-fluid-penetrated soft- 
particles 0 as a function of the Debye length for the 


two applied fields are both included in the figure. For 
weak fields, there is a stronger dependence on Debye 
length, and the electrophoretic mobility approaches to 
zero, meaning that the fluctuating hydrodynamic inter¬ 
actions (driven by Brownian motion) affect particle tra¬ 
jectories, averaging down the electrophoretic induced mo¬ 
tion. On the other hand, for strong fields, the normalized 
mobility collapses into a single curve, maintaining the 
dominance of the electrokinetic forces. Ohshima’s solu¬ 
tion is obtained from the linearized PB; other authors 
had suggested different solutions that result in similar 
approximations for these electrokinetic properties, which 
are restricted to zero concentrations and weak potentials. 
The approach we proposed here avoids the approximation 
of electrokinetic properties and relationships. It provides 
a platform for other colloidal and polymeric systems, in¬ 
cluding deformable objects, rigid suspensions and cells, 
where lubrication forces and slip conditions are impor¬ 
tant. 

With the method described here, it is possible to per¬ 
form efficient simulations of discrete charged elements 
embedded in charged solvents at non-equilibrium con¬ 
ditions. It should And applications in a wide variety of 
situations involving the physics of polymeric and colloidal 
systems including micro- and nano-scale systems. 
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